set more off	
	
set seed 06031974

	forval b= 1(1)100 {

	di `b'
	
	use Data/Andrea/PV_forboot.dta, clear
	
	bsample , cluster(sid)
	
	
	quietly reg perkW_energy_generate i(6/18).hour##c.temp_interp i(6/18).hour##c.ghi i(6/18).hour#i(2/4).season##c.latitude i(6/18).hour#i(2/4).season##c.longitude
	eststo reg4f

use Data/Andrea/sample_forboot_PV.dta, clear


	bsample , cluster(account_number) idcluster(id)
	drop household_hour
	quietly predict PV_predictionkWh if hour >= 6 & hour<= 19 


** censor PV production at zero
	quietly gen PV_predict_cenkWh = PV_prediction 
	quietly replace PV_predict_cenkWh = 0 if PV_prediction <0

	quietly gen PV_productionkWh = PV_predict_cen* sizekW

	quietly gen unobserved_use_kWh = PV_production - exportkWh

	quietly gen consumption_kWh =  unobserved_use_kWh + importkWh
	quietly replace consumption = importkWh if unobserved_use_kWh ==.

	quietly replace PV_production = 0 if PV_production ==. &( hour<6 | hour>19)

	quietly gen virtualincome = price * PV_production

	keep PV_production consumption_kWh virtualincome exportkWh pexport_FIT price hours* weekend hour
	
	gen price_coef =.
	gen export_elasticity=.
	gen consumption_elasticity =.
	gen PVlevel=.

	save Data/current_bsample.dta, replace
	
	
	foreach model in Baselineprice  {
	
	use Data/current_bsample.dta, clear
	estimates use Results/Andrea/`model'_`b'

 
	replace price_coef= .
 
	forval x = 1(1)24 {
	replace price_coef = _b[hour_price`x'] if hours`x' ==1
 
					}
	
	
		
	
		}
	
	
	replace export_elasticity = (pexport_FIT/exportkWh)* (-price_coef ) if exportkWh >0.001 & hour >=6 & hour< = 18
	replace consumption_elasticity = (price/consumption_kWh)* (price_coef ) if consumption_kWh >0.001
	
	
	** collapse by hour by PV for plotting (three levels of PV)? or maybe just low and high as in Figure 6.  low < 0.5, med = 0.5-3, high >= 3

	replace PVlevel = (PV_production < 0.5)
	replace PVlevel = 2 if PV_production >=0.5 & PV_production< 3.5
	replace PVlevel = 3 if PV_production >=3.5 & PV_production<. 
	
	collapse (mean) consumption_elasticity export_elasticity , by(hour PVlevel)
	
	if `b'==1 {
	
	save Results/Andrea/`model'_elasticity_noinc, replace
	
			}
	
	else {
	
	append using Results/Andrea/`model'_elasticity_noinc
	
	save Results/Andrea/`model'_elasticity_noinc, replace
	
			}
	
	
													}
	
				
use Results/Andrea/_elasticity_noinc.dta, clear
	
	
** collapse and make an area graph when have all the results
	collapse (mean) consumption_elasticity export_elasticity (sd) ///
	consumption_elasticity_sd = consumption_elasticity export_elasticity_sd = export_elasticity, by(PVlevel hour)

	foreach var in consumption_elasticity export_elasticity {
	gen `var'_high = `var' + 1.96*`var'_sd
	gen `var'_low  = `var' - 1.96*`var'_sd
	}


	
	twoway (rarea consumption_elasticity_high consumption_elasticity_low hour if PVlevel==1, color(gs10)  lcolor(gs10)) ///
	(rarea consumption_elasticity_high consumption_elasticity_low hour if PVlevel==2 , color(bluishgray%80))  ///
	(rarea consumption_elasticity_high consumption_elasticity_low hour if PVlevel==3, color(eltgreen%60) )  /// 
	(connected consumption_elasticity hour if PVlevel==1, lcolor(gs8) msymbol(none) lpattern(dash) lwidth(thin)) ///
	(connected consumption_elasticity hour if PVlevel==2 , lcolor(edkblue)  msymbol(none) lpattern(shortdash) lwidth(thin) )  ///
	(connected consumption_elasticity hour if PVlevel==3, lcolor(teal)  msymbol(none) lpattern(solid) lwidth(thin) lwidth(thin)),  /// 
	legend(order(4 "Low production" 5 "Medium production" 6 "High production") region(style(none) lcolor(none))) ///
	 ytitle("Consumption elasticity")  graphregion(color(white)) ///
	xscale(range(0 23))  xlabel(0(4)20)
graph export Results/Andrea/noinc_consumptionelasticity.pdf, replace




	twoway (rarea export_elasticity_high export_elasticity_low hour if PVlevel==1, color(gs10) lcolor(gs10) ) ///
	(rarea export_elasticity_high export_elasticity_low hour if PVlevel==2 , color(bluishgray%80) lcolor(bluishgray))  ///
	(rarea export_elasticity_high export_elasticity_low hour if PVlevel==3, color(eltgreen%60) lcolor(eltgreen))  /// 
	(connected export_elasticity hour if PVlevel==1, lcolor(gs8) msymbol(none)  lpattern(dash) lwidth(thin) ) ///
	(connected export_elasticity hour if PVlevel==2 , lcolor(edkblue)  msymbol(none) lpattern(shortdash) lwidth(thin) )  ///
	(connected export_elasticity hour if PVlevel==3, lcolor(teal)  msymbol(none) lpattern(solid) lwidth(thin)) if hour>=6 & hour<=18,  /// 
	legend(order(4 "Low production" 5 "Medium production" 6 "High production") region(style(none) lcolor(none))) ///
	 ytitle("Export elasticity")  graphregion(color(white)) ///
	xscale(range(6 18)) xlabel(6(3)18)
graph export Results/Andrea/noinc_exportelasticity.pdf, replace

													

	
	
	
